PETupdate Subroutine

SUBROUTINE PETupdate &
!
( time )

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT(IN) :: time

!local declarations:
INTEGER (KIND = short)   :: i, j
INTEGER (KIND = short)   :: jDay !!julian day
REAL (KIND = float)      :: sunsetHourAngle !!sunset hour angle [radians]
REAL (KIND = float)      :: relDist !!relative distance between the earth and the sun
REAL (KIND = float)      :: declination !! declination
REAL (KIND = float)      :: etRad !!extra terrestrial solar radiation [mm/day]


!------------------------------end of declarations-----------------------------

!Update Kc
IF (dtKc > 0 ) THEN
	IF( time == tNewKc ) THEN
		CALL  KcUpdate (time)
		tNewKc = tNewKc + dtKc
	END IF
END IF	


IF ( time == tnewET )  THEN !it's time to update ET
    
    !set next time to update PET
    tnewET = tnewET + dtET
    
    IF ( compute_hargreaves ) THEN !compute values for applying Hargreaves
        
        !calculate solar declination
        jDay = DayOfYear(time,'noleap') 
        declination = 0.4093 * SIN( 2. * pi * jDay / 365. - 1.405 )
        
        !calculate sunset hour angle
        sunsetHourAngle = ACOS( - TAN (latCentroid) * TAN (declination) )
        
        !calculate relative distance between the earth and the sun
        relDist = 1. + 0.033 * COS(2. * pi * jDay / 365.)
        
        !calculate extra terrestrial solar radiation (Allen et al. 1998)
        etRad = 15.392 * relDist * ( sunsetHourAngle * SIN(latCentroid) * &
                                     SIN(declination) + COS(latCentroid) * &
                                    COS(declination) * SIN(sunsetHourAngle) )
        
    END IF
    
   
    
    DO i = 1, model_map % idim
        DO j = 1, model_map % jdim
            IF ( model_map % mat (i,j) /= model_map % nodata) THEN
                SELECT CASE (model_map % mat (i,j) )
                    
                CASE (HARGREAVES)
                   CALL HargreavesSamani (time, temperatureAirDailyMean % mat(i,j), &
                                         temperatureAirDailyMax % mat(i,j), &
                                         temperatureAirDailyMin % mat (i,j), &
                                         etRad, pet % mat(i,j) )
                   

                   pe % mat (i,j) = pet % mat(i,j)
                   
                   pt % mat (i,j) = pet % mat(i,j)
                   
                CASE (HARGREAVES_MOD)
                   CALL HargreavesSamaniModified (time, &
                                         temperatureAirDailyMean % mat(i,j), &
                                         temperatureAirDailyMax % mat(i,j), &
                                         temperatureAirDailyMin % mat (i,j), &
                                         etRad, dem % mat (i,j), &
                                         pet % mat(i,j) )
                   

                   pe % mat (i,j) = pet % mat(i,j)
                   
                   pt % mat (i,j) = pet % mat(i,j)
                   
                CASE (PRIESTLEY_TAYLOR)  
                    
                   CALL PriestleyTaylor (temperatureAir % mat (i,j), &
                                         netRadiation % mat (i,j), &
                                         fvcover % mat (i,j), &
                                         dem % mat (i,j), &
                                         pt % mat (i,j), &
                                         pe % mat (i,j), &
                                         pet % mat (i,j) )
                   
                CASE (PENMAN_MONTEITH) 
                    
                    
                    
               !CALL  PetPM (R = netRadiation % mat (i,j), &
               !            T = temperatureAir % mat (i,j), &
               !            W =  windSpeed % mat (i,j), &
               !            rs = rsMin % mat (i,j), &
               !            v = fvcover % mat (i,j), &
               !            lai = lai % mat (i,j), &
               !            z = plantsHeight % mat (i,j), &
               !            um = relHumidityAir % mat (i,j), &
               !            quota = dem % mat (i,j), &
               !            pet = pet % mat (i,j), &
               !            pe = pe % mat (i,j), &
               !            pt = pt % mat (i,j))
                    
                    
                    
                   CALL PenmanMonteith (temperatureAir % mat (i,j), &
                                         netRadiation % mat (i,j), &
                                         relHumidityAir % mat (i,j), &
                                         windSpeed % mat (i,j), &
                                         fvcover % mat (i,j), &
                                         plantsHeight % mat (i,j), &
                                         dem % mat (i,j), &
                                         anemometersWS % offsetZ, &
                                         hygrometers % offsetZ, &
                                         lai % mat (i,j), &
                                         rsMin % mat (i,j), &
                                         pt % mat (i,j), &
                                         pe % mat (i,j), &
                                         pet % mat (i,j) )
                   
                 CASE (FAO56_PENMAN_MONTEITH)  
                    
                   CALL FAO56PenmanMonteith (temperatureAir % mat (i,j), &
                                         netRadiation % mat (i,j), &
                                         netRadiationFAO % mat (i,j), &
                                         relHumidityAir % mat (i,j), &
                                         windSpeed % mat (i,j), &
                                         fvcover % mat (i,j), &
                                         dem % mat (i,j), &
                                         anemometersWS % offsetZ, &
                                         hygrometers % offsetZ, &
                                         lai % mat (i,j), &
                                         pt % mat (i,j), &
                                         pe % mat (i,j), &
                                         pet % mat (i,j) )
  
                 END SELECT
                 
                 
                 IF ( useCropCoefficient) THEN
                       pt % mat (i,j) = pet % mat(i,j) * kc_map % mat (i,j)
                       !update pet
                       pet % mat (i,j) = fvcover % mat (i,j) * pt % mat (i,j) + &
                                        ( 1. - fvcover % mat (i,j) ) *  pe % mat (i,j) 
                 ELSE
                       pt % mat (i,j) = pet % mat(i,j)
                 END IF
                 
                 !check negative values
                 IF ( pt % mat (i,j) < 0.) pt % mat (i,j) = 0.
                 IF ( pe % mat (i,j) < 0.) pe % mat (i,j) = 0.
                 IF ( pet % mat (i,j) < 0.) pet % mat (i,j) = 0.
                 
                 
            END IF
        END DO
    END DO
    
END IF

RETURN
END SUBROUTINE PETupdate